function svl_jpe_figures()
% Sections 3 & 4, Appendix B.2

clear
clc
close all

[ell,theta,sigma,pi,rho,eta,i,M,A,B]=parameters();

% figure2(ell,theta,sigma,pi,rho,eta,i,M,A,B)    %   6.046038 seconds
% figure3(ell,theta,sigma,pi,rho,eta,i,M,A,B)    %   4.795128 seconds
% figure4(ell,theta,sigma,pi,rho,eta,i,M,A,B)    %  27.073627 seconds
% figure5(ell,theta,sigma,pi,rho,eta,i,M,A,B)    % 308.779134 seconds
% figure6(ell,theta,sigma,pi,rho,eta,i,M,A,B)    % 323.064565 seconds
% figure7(ell,theta,sigma,pi,rho,eta,i,M,A,B)    %  51.461650 seconds
% figure8(ell,theta,sigma,pi,rho,eta,i,M,A,B)    % 297.998642 seconds
% figure10(ell,theta,sigma,pi,rho,eta,i,M,A,B)   % 859.837159 seconds
% figure11(ell,theta,sigma,pi,rho,eta,i,M,A,B)   %  27.315433 seconds

%=========================================================================
% parameters
%-------------------------------------------------------------------------
function [ell,theta,sigma,pi,rho,eta,i,M,A,B]=parameters()
% parameter setting
ell   = 0.5;       % probability of being C-type
theta = 0.5;       % C-type's bargaining power
sigma = 1;         % CRRA IES
pi    = 0.95;      % 1 - default rate
rho   = 0;         % recovery rate
eta   = 0;         % IRS
i     = 0.1;
M     = 1;
A     = 0.05;
B     = 0.05;

%=========================================================================
% subfunctions
%-------------------------------------------------------------------------
function u = u(q,sigma)
% CRRA utility: u(q)
if sigma < 1
    u = q.^(1-sigma)./(1-sigma);
elseif sigma == 1
    u = log(q);
end

function mu = mu(q,sigma)
% marginal utility: u'(q)
mu = q.^(-sigma);

function w = w(q,sigma,theta)
% w function in effective bargaining power
w = theta + (1-theta)*mu(q,sigma);

%=========================================================================
% solvers
%-------------------------------------------------------------------------
function [q0A,q1A,q0B,q1Bn,eNn]=symmetricEquilibrium(eC,ell,theta,sigma,pi,rho,eta,i,M,A,B)
% partial equilibrium, given eC

lb = [0,0,0,0,0.5]; % lower bound constraint
ub = [1,1,1,1,0.5]; % upper bound constraint
x0 = [0.8,0.9,0.8,0.9,0.5]; % initial guess
opts = optimoptions(@fmincon,'Algorithm','sqp','Display','off'); 
nonlcon  = @(x) fminconstr(x,eC,ell,theta,sigma,pi,rho,eta,i,M,A,B);
x = fmincon(@(x)0,x0,[],[],[],[],lb,ub,nonlcon,opts);

q0A  = x(1);
q1A  = x(2);
q0B  = x(3);
q1Bn = x(4);
eNn  = x(5);
eNd  = 1;

function [G_eC,L_A,L_B,q0A,q1A,q0B,q1Bn,eNn]=partialEquilibrium(eC,ell,theta,sigma,pi,rho,eta,i,M,A,B)
% partial equilibrium, given eC

lb = [0,0,0,0,0]; % lower bound constraint
ub = [1,1,1,1,1]; % upper bound constraint
x0 = [0.8,0.9,0.8,0.9,A/(A+B)]; % initial guess
opts = optimoptions(@fmincon,'Algorithm','sqp','Display','off'); 
if eC > 0 && eC < 1
    nonlcon  = @(x) fminconstr(x,eC,ell,theta,sigma,pi,rho,eta,i,M,A,B);
    x = fmincon(@(x)0,x0,[],[],[],[],lb,ub,nonlcon,opts);
elseif eC == 0
    nonlcon0 = @(x) fminconstr0(x,eC,ell,theta,sigma,pi,rho,eta,i,M,A,B);
    x = fmincon(@(x)0,x0,[],[],[],[],lb,ub,nonlcon0,opts);
elseif eC == 1
    nonlcon1 = @(x) fminconstr1(x,eC,ell,theta,sigma,pi,rho,eta,i,M,A,B);
    x = fmincon(@(x)0,x0,[],[],[],[],lb,ub,nonlcon1,opts);
end

q0A  = x(1);
q1A  = x(2);
q0B  = x(3);
q1Bn = x(4);
eNn  = x(5);
eNd  = 1;

if eC > 0 && eC < 1
    a_CA_n =     eNn*(1-ell)/(eC*ell+eNn*(1-ell))^(1-eta);
    a_CB_n = (1-eNn)*(1-ell)/((1-eC)*ell+(1-eNn)*(1-ell))^(1-eta);
    a_CA_d =     eNd*(1-ell)/(eC*ell+eNd*(1-ell))^(1-eta);
	S_CA   = theta*(u(q1A,sigma) -u(q0A,sigma)-q1A +q0A);
	S_CB_n = theta*(u(q1Bn,sigma)-u(q0B,sigma)-q1Bn+q0B);
	L_A         = ell*(pi*a_CA_n+(1-pi)*a_CA_d)*theta/w(q1A,sigma,theta)*(mu(q1A,sigma)-1);
	L_B_rho_bar = ell*pi*a_CB_n*theta/w(q1Bn,sigma,theta)*(mu(q1Bn,sigma)-1);
	S_CA_tilde  = - i*q0A - L_A*((1-theta)*(u(q1A,sigma)-u(q0A,sigma))+theta*(q1A-q0A)) + ell*(u(q0A,sigma)-q0A+(pi*a_CA_n+(1-pi)*a_CA_d)*S_CA);
	S_CB_tilde  = - i*q0B - L_B_rho_bar*((1-theta)*(u(q1Bn,sigma)-u(q0B,sigma))+theta*(q1Bn-q0B)) + ell*(u(q0B,sigma)-q0B+pi*a_CB_n*S_CB_n);
	G_eC = S_CA_tilde - S_CB_tilde;
    L_A  = 100 * L_A;
    L_B  = 100 * L_B_rho_bar / pi;
elseif eC == 0
	a_CA_n = 0;
	a_CB_n = 1-ell;
	a_CA_d = 0;
	S_CA   = theta*(u(q1A,sigma) -u(q0A,sigma)-q1A +q0A);
	S_CB_n = theta*(u(q1Bn,sigma)-u(q0B,sigma)-q1Bn+q0B);
	L_A         = 0;
	L_B_rho_bar = ell*pi*a_CB_n*theta/w(q1Bn,sigma,theta)*(mu(q1Bn,sigma)-1);
	S_CA_tilde  = - i*q0A + ell*(u(q0A,sigma)-q0A);
	S_CB_tilde  = - i*q0B - L_B_rho_bar*((1-theta)*(u(q1Bn,sigma)-u(q0B,sigma))+theta*(q1Bn-q0B)) + ell*(u(q0B,sigma)-q0B+pi*a_CB_n*S_CB_n);
	G_eC = S_CA_tilde - S_CB_tilde;
    L_A  = 100 * L_A;
    L_B  = 100 * L_B_rho_bar / pi;
elseif eC == 1
	a_CA_n = 1-ell;
	a_CB_n = 0;
	a_CA_d = 1-ell;
	S_CA   = theta*(u(q1A,sigma) -u(q0A,sigma)-q1A +q0A);
	S_CB_n = theta*(u(q1Bn,sigma)-u(q0B,sigma)-q1Bn+q0B);
	L_A         = ell*(pi*a_CA_n+(1-pi)*a_CA_d)*theta/w(q1A,sigma,theta)*(mu(q1A,sigma)-1);
	L_B_rho_bar = 0;
	S_CA_tilde  = - i*q0A - L_A*((1-theta)*(u(q1A,sigma)-u(q0A,sigma))+theta*(q1A-q0A)) + ell*(u(q0A,sigma)-q0A+(pi*a_CA_n+(1-pi)*a_CA_d)*S_CA);
	S_CB_tilde  = - i*q0B + ell*(u(q0B,sigma)-q0B);
	G_eC = S_CA_tilde - S_CB_tilde;
    L_A  = 100 * L_A;
    L_B  = 100 * L_B_rho_bar / pi;
end

function [eC,L_A,L_B,q0A,q1A,q0B,q1Bn,eNn]=fixedPoint(ell,theta,sigma,pi,rho,eta,i,M,A,B)
% general equilibrium

% initial guess
eC = A/(A+B);

cnt = 0;
precision = 1;
while precision > 10^-6
    cnt = cnt + 1;
    if cnt > 1
        eC = temp_eC;
    end
    
    fprintf('finding a fixed point ... iter: %.0f ... ell = %.2f, theta = %.2f, eta = %.2f, i = %.4f, pi = %.5f, A= %.4f, B = %.4f ... eC = %.4f ',cnt,ell,theta,eta,i,pi,A,B,eC)
    [G_eC,L_A,L_B,q0A,q1A,q0B,q1Bn,eNn]=partialEquilibrium(eC,ell,theta,sigma,pi,rho,eta,i,M,A,B);
    fprintf('... G_eC = %.4f\n\n',G_eC)
    
    itr = 200; %1
    if cnt < itr
        if eC < 1
            if eC > 0
%                 temp_eC = max(0,min(1,eC + .5*G_eC)); 
                temp_eC = max(0,min(1,eC + 5*G_eC)); 
            elseif eC == 0
                break
            end
        elseif eC == 1
            break
        end
        precision = abs(G_eC);
    elseif cnt >= itr
        if cnt > itr
            temp_eC_2 = temp_eC_1;
        end
        if eC < 1
            if eC > 0
                temp_eC_1 = eC;
                temp_eC = max(0,min(1,eC + sign(G_eC)*10^-4)); 
            elseif eC == 0
                break
            end
        elseif eC == 1
            break
        end
        if cnt > itr
            if temp_eC_2 == temp_eC
                break
            end
        end
    end
end

%=========================================================================
% grid
%-------------------------------------------------------------------------
function [GeC,EC,LA,LB,Q0A,Q1A,Q0B,Q1Bn,ENn]=EC_mat(ell,theta,sigma,pi,rho,eta,i,M,A,B)
% partial equilibrium, given eC

% EC = 0:0.01:1;
EC = [0 0.001 0.01:0.01:0.99 0.999 1]; 

GeC  = zeros(1,length(EC));
LA   = zeros(1,length(EC));
LB   = zeros(1,length(EC));
Q0A  = zeros(1,length(EC));
Q1A  = zeros(1,length(EC));
Q0B  = zeros(1,length(EC));
Q1Bn = zeros(1,length(EC));
ENn  = zeros(1,length(EC));

for j=1:length(EC)
    eC = EC(j);
    
    fprintf('Computing ... eC = %.4f\n\n',eC)
    [G_eC,L_A,L_B,q0A,q1A,q0B,q1Bn,eNn]=partialEquilibrium(eC,ell,theta,sigma,pi,rho,eta,i,M,A,B);
    
    GeC(j)  = G_eC;
    LA(j)   = L_A;
    LB(j)   = L_B;
    Q0A(j)  = q0A;
    Q1A(j)  = q1A;
    Q0B(j)  = q0B;
    Q1Bn(j) = q1Bn;
    ENn(j)  = eNn;
end

function [PI,EC,LA,LB,Q0A,Q1A,Q0B,Q1Bn,ENn]=PI_mat(ell,theta,sigma,pi,rho,eta,i,M,A,B)
% general equilibrium: change pi, given A,B

PI = 0.9:0.01:1; 

EC   = zeros(1,length(PI));
LA   = zeros(1,length(PI));
LB   = zeros(1,length(PI));
Q0A  = zeros(1,length(PI));
Q1A  = zeros(1,length(PI));
Q0B  = zeros(1,length(PI));
Q1Bn = zeros(1,length(PI));
ENn  = zeros(1,length(PI));

for j=1:length(PI)
    pi = PI(j);
    
    [eC,L_A,L_B,q0A,q1A,q0B,q1Bn,eNn]=fixedPoint(ell,theta,sigma,pi,rho,eta,i,M,A,B);
    
    EC(j)   = eC;
    LA(j)   = L_A;
    LB(j)   = L_B;
    Q0A(j)  = q0A;
    Q1A(j)  = q1A;
    Q0B(j)  = q0B;
    Q1Bn(j) = q1Bn;
    ENn(j)  = eNn;
end

function [Bmat,EC,LA,LB,Q0A,Q1A,Q0B,Q1Bn,ENn]=B_mat(ell,theta,sigma,pi,rho,eta,i,M,A,B)
% general equilibrium: change B, given A

global Bmin Bgap Bmax 
Bmat = Bmin:Bgap:Bmax; 

EC   = zeros(1,length(Bmat));
LA   = zeros(1,length(Bmat));
LB   = zeros(1,length(Bmat));
Q0A  = zeros(1,length(Bmat));
Q1A  = zeros(1,length(Bmat));
Q0B  = zeros(1,length(Bmat));
Q1Bn = zeros(1,length(Bmat));
ENn  = zeros(1,length(Bmat));

for j=1:length(Bmat)
    B = Bmat(j);
    
    [eC,L_A,L_B,q0A,q1A,q0B,q1Bn,eNn]=fixedPoint(ell,theta,sigma,pi,rho,eta,i,M,A,B);
    
    EC(j)   = eC;
    LA(j)   = L_A;
    LB(j)   = L_B;
    Q0A(j)  = q0A;
    Q1A(j)  = q1A;
    Q0B(j)  = q0B;
    Q1Bn(j) = q1Bn;
    ENn(j)  = eNn;
end

function [Amat,EC,LA,LB,Q0A,Q1A,Q0B,Q1Bn,ENn]=A_mat(ell,theta,sigma,pi,rho,eta,i,M,A,B)
% general equilibrium: change A, given B

global Amin Agap Amax 
Amat = [Amin:Agap:Amax Amax]; 

EC   = zeros(1,length(Amat));
LA   = zeros(1,length(Amat));
LB   = zeros(1,length(Amat));
Q0A  = zeros(1,length(Amat));
Q1A  = zeros(1,length(Amat));
Q0B  = zeros(1,length(Amat));
Q1Bn = zeros(1,length(Amat));
ENn  = zeros(1,length(Amat));

for j=1:length(Amat)
    A = Amat(j);
    
    [eC,L_A,L_B,q0A,q1A,q0B,q1Bn,eNn]=fixedPoint(ell,theta,sigma,pi,rho,eta,i,M,A,B);
    
    EC(j)   = eC;
    LA(j)   = L_A;
    LB(j)   = L_B;
    Q0A(j)  = q0A;
    Q1A(j)  = q1A;
    Q0B(j)  = q0B;
    Q1Bn(j) = q1Bn;
    ENn(j)  = eNn;
end

%=========================================================================
% subplots
%-------------------------------------------------------------------------
function p=sub_figure2_GeC(ell,theta,sigma,pi,rho,eta,i,M,A,B,color)

[GeC,EC,LA,LB,Q0A,Q1A,Q0B,Q1Bn,ENn]=EC_mat(ell,theta,sigma,pi,rho,eta,i,M,A,B);

p=plot(EC(2:end-1),GeC(2:end-1),'Color',color,'Linewidth',3); hold on
plot(EC(1),GeC(1),'o','Color',color,'Linewidth',3), hold on
plot(EC(end),GeC(end),'o','Color',color,'Linewidth',3), hold on
plot(EC,zeros(1,length(EC))','k--','Linewidth',2), hold on

function p=sub_figure2_entry(ell,theta,sigma,pi,rho,eta,i,M,A,B,color)

[GeC,EC,LA,LB,Q0A,Q1A,Q0B,Q1Bn,ENn]=EC_mat(ell,theta,sigma,pi,rho,eta,i,M,A,B);

p=plot(EC(2:end-1),ENn(2:end-1),'Color',color,'Linewidth',3); hold on
plot(EC(1),ENn(1),'o','Color',color,'Linewidth',3), hold on
plot(EC(end),ENn(end),'o','Color',color,'Linewidth',3), hold on

function [PI,DELTA_A,DELTA_B]=sub_figure11(ell,theta,sigma,pi,rho,eta,i,M,A,B)

[PI,EC,LA,LB,Q0A,Q1A,Q0B,Q1Bn,ENn]=PI_mat(ell,theta,sigma,pi,rho,eta,i,M,A,B);

eNd = 1;
aCAn =     ENn*(1-ell)./(EC*ell+ENn*(1-ell)).^(1-eta);
aCBn = (1-ENn)*(1-ell)./((1-EC)*ell+(1-ENn)*(1-ell)).^(1-eta);
aCAd =     eNd*(1-ell)./(EC*ell+eNd*(1-ell)).^(1-eta);

DELTA_A_N = EC*ell.*aCAn.*min(M*((1-theta)*(log(1)-log(Q0A))+theta*(1-Q0A))./(EC.*Q0A+(1-EC).*Q0B),A./EC);
DELTA_A_D = EC*ell.*aCAd.*min(M*((1-theta)*(log(1)-log(Q0A))+theta*(1-Q0A))./(EC.*Q0A+(1-EC).*Q0B),A./EC);
DELTA_A   = PI.*DELTA_A_N+(1-PI).*DELTA_A_D;

DELTA_B_N = (1-EC)*ell.*aCBn.*min(M*((1-theta)*(log(1)-log(Q0B))+theta*(1-Q0B))./(EC.*Q0A+(1-EC).*Q0B),B./(1-EC));
DELTA_B   = PI.*DELTA_B_N;

function sub_figure5(ell,theta,sigma,pi,rho,eta,i,M,A,B)

[Bmat,EC,LA,LB,Q0A,Q1A,Q0B,Q1Bn,ENn]=B_mat(ell,theta,sigma,pi,rho,eta,i,M,A,B);

plot(Bmat,LA,'k-','Linewidth',3), hold on
plot(Bmat,LB,'k--','Linewidth',3), hold on, 
plot([A,A],[0,4],'k--','Linewidth',2,'Color',[0.5 0.5 0.5])
axis square
set(gca,'FontSize',12,'TickLabelInterpreter','latex')
ylim([0 1.8])
xticks([0 0.05 0.1 0.2 0.3])
xticklabels({'0','$$S_A$$','0.1','0.2','0.3'})

function sub_figure6(ell,theta,sigma,pi,rho,eta,i,M,A,B)

[Bmat,EC,LA,LB,Q0A,Q1A,Q0B,Q1Bn,ENn]=B_mat(ell,theta,sigma,pi,rho,eta,i,M,A,B);

aCBn = (1-ENn)*(1-ell)./((1-EC)*ell+(1-ENn)*(1-ell)).^(1-eta);
aNBn =      (1-EC)*ell./((1-EC)*ell+(1-ENn)*(1-ell)).^(1-eta);

plot(Bmat,aCBn,'k-','Linewidth',3), hold on
plot(Bmat,aNBn,'k--','Linewidth',3), hold on
plot([A,A],[0.175 0.525],'--','Linewidth',2,'Color',[0.5 0.5 0.5])
axis square
set(gca,'FontSize',12,'TickLabelInterpreter','latex')
ylim([0.175 0.525])
xticks([0 0.05 0.1 0.2 0.3])
xticklabels({'0','$$S_A$$','0.1','0.2','0.3'})

function sub_figure8(ell,theta,sigma,pi,rho,eta,i,M,A,B)

[Amat,EC,LA,LB,Q0A,Q1A,Q0B,Q1Bn,ENn]=A_mat(ell,theta,sigma,pi,rho,eta,i,M,A,B);

aCAn = zeros(1,length(EC));
aCBn = zeros(1,length(EC));
aCAd = zeros(1,length(EC));
for j=1:length(EC)
    eC  = EC(j);
    eNn = ENn(j);
    eNd = 1;
    if eC > 0 && eC < 1
        aCAn(j) =     eNn*(1-ell)/(eC*ell+eNn*(1-ell))^(1-eta);
        aCBn(j) = (1-eNn)*(1-ell)/((1-eC)*ell+(1-eNn)*(1-ell))^(1-eta);
        aCAd(j) =     eNd*(1-ell)/(eC*ell+eNd*(1-ell))^(1-eta);
    elseif eC == 0
        aCAn(j) = 0;
        aCBn(j) = 1-ell;
        aCAd(j) = 0;
    elseif eC == 1
        aCAn(j) = 1-ell;
        aCBn(j) = 0;
        aCAd(j) = 1-ell;
    end
end

WELFARE_N  = EC*ell.*((1-aCAn).*(log(Q0A)-Q0A)+aCAn.*(log(Q1A)-Q1A))+(1-EC)*ell.*((1-aCBn).*(log(Q0B)-Q0B)+aCBn.*(log(Q1Bn)-Q1Bn));
WELFARE_D  = EC*ell.*((1-aCAd).*(log(Q0A)-Q0A)+aCAd.*(log(Q1A)-Q1A))+(1-EC)*ell.*(log(Q0B)-Q0B);
WELFARE    = pi*WELFARE_N+(1-pi)*WELFARE_D;

plot(Amat,WELFARE,'k-','Linewidth',3), hold on
plot([B,B],[-0.5071 -0.5059],'k--','Linewidth',2,'Color',[0.5 0.5 0.5])
axis square
set(gca,'FontSize',12,'TickLabelInterpreter','latex')
ylim([-0.5071 -0.5059]), yticks([-0.507 -0.5068 -0.5066 -0.5064 -0.5062 -0.506]), yticklabels({' ',' ',' ',' ',' ',' '})
xticks([0 0.05 0.1 0.2 0.3]), xticklabels({'0','$$S_B$$','0.1','0.2','0.3'})

%=========================================================================
% figures 2, 3
%-------------------------------------------------------------------------
function figure2(ell,theta,sigma,pi,rho,eta,i,M,A,B)

eta = 0;

figure()
cmap=colormap(gray);
c1 = cmap(40,:);
c2 = cmap(25,:);
c3 = cmap(1,:);

subaxis(1,2,1, 'Spacing', 0.06, 'Padding', 0, 'Margin', 0.08);
% small A
A = .05; B = .14;
p1=sub_figure2_GeC(ell,theta,sigma,pi,rho,eta,i,M,A,B,c1);
p1.LineStyle = ':';
% medium A
A = .09; B = .14;
p2=sub_figure2_GeC(ell,theta,sigma,pi,rho,eta,i,M,A,B,c2);
p2.LineStyle = '--';
% large A
A = .14; B = .14;
p3=sub_figure2_GeC(ell,theta,sigma,pi,rho,eta,i,M,A,B,c3);
p3.LineStyle = '-';

axis square
xticks(0:0.2:1)
ylim([-0.006 0.006]), yticks([-0.006 -0.004 -0.002 0 0.002 0.004 0.006]), yticklabels({' ',' ',' ','0',' ',' ',' '})
set(legend([p3(1) p2(1) p1(1)],{'\,large $$S_A$$','\,medium $$S_A$$','\,small $$S_A$$'},'Interpreter','latex'), 'FontSize', 13,'Location','southwest')
set(gca,'FontSize',10.5,'TickLabelInterpreter','latex')
set(xlabel('$$e_C^{}$$','interpreter','latex'), 'FontSize', 16)
set(title('$$G(e_C^{})$$','interpreter','latex'), 'FontSize', 16)

subaxis(1,2,2, 'Spacing', 0.06, 'Padding', 0, 'Margin', 0.08);
% small A
A = .05; B = .14;
p1=sub_figure2_entry(ell,theta,sigma,pi,rho,eta,i,M,A,B,c1);
p1.LineStyle = ':';
% medium A
A = .09; B = .14;
p2=sub_figure2_entry(ell,theta,sigma,pi,rho,eta,i,M,A,B,c2);
p2.LineStyle = '--';
% large A
A = .14; B = .14;
p3=sub_figure2_entry(ell,theta,sigma,pi,rho,eta,i,M,A,B,c3);
p3.LineStyle = '-';

axis square
xticks(0:0.2:1), yticks(0:0.2:1)
set(gca,'FontSize',10.5,'TickLabelInterpreter','latex')
set(xlabel('$$e_C^{}$$','interpreter','latex'), 'FontSize', 16)
set(title('$$e_N^n$$','interpreter','latex'), 'FontSize', 16)

set(gcf,'pos',[30 70 830 430])
% print('fig2','-depsc')

disp('done')

function figure3(ell,theta,sigma,pi,rho,eta,i,M,A,B)

eta = 0.3;

figure()
cmap=colormap(gray);
c1 = cmap(40,:);
c2 = cmap(25,:);
c3 = cmap(1,:);

subaxis(1,2,1, 'Spacing', 0.06, 'Padding', 0, 'Margin', 0.08);
% small A
A = .05; B = .05;
p1=sub_figure2_GeC(ell,theta,sigma,pi,rho,eta,i,M,A,B,c1);
p1.LineStyle = ':';
% medium A
A = .09; B = .05;
p2=sub_figure2_GeC(ell,theta,sigma,pi,rho,eta,i,M,A,B,c2);
p2.LineStyle = '--';
% large A
A = .14; B = .05;
p3=sub_figure2_GeC(ell,theta,sigma,pi,rho,eta,i,M,A,B,c3);
p3.LineStyle = '-';

axis square
xticks(0:0.2:1)
ylim([-0.003 0.004]), yticks([-0.003 -0.002 -0.001 0 0.001 0.002 0.003 0.004]), yticklabels({' ',' ',' ','0',' ',' ',' ',' '})
set(legend([p3(1) p2(1) p1(1)],{'\,large $$S_A$$','\,medium $$S_A$$','\,small $$S_A$$'},'Interpreter','latex'), 'FontSize', 13,'Location','southwest')
set(gca,'FontSize',10.5,'TickLabelInterpreter','latex')
set(xlabel('$$e_C^{}$$','interpreter','latex'), 'FontSize', 16)
set(title('$$G(e_C^{})$$','interpreter','latex'), 'FontSize', 16)

subaxis(1,2,2, 'Spacing', 0.06, 'Padding', 0, 'Margin', 0.08);
% small A
A = .05; B = .05;
p1=sub_figure2_entry(ell,theta,sigma,pi,rho,eta,i,M,A,B,c1);
p1.LineStyle = ':';
% medium A
A = .09; B = .05;
p2=sub_figure2_entry(ell,theta,sigma,pi,rho,eta,i,M,A,B,c2);
p2.LineStyle = '--';
% large A
A = .14; B = .05;
p3=sub_figure2_entry(ell,theta,sigma,pi,rho,eta,i,M,A,B,c3);
p3.LineStyle = '-';

axis square
xticks(0:0.2:1), yticks(0:0.2:1)
set(gca,'FontSize',10.5,'TickLabelInterpreter','latex')
set(xlabel('$$e_C^{}$$','interpreter','latex'), 'FontSize', 16)
set(title('$$e_N^n$$','interpreter','latex'), 'FontSize', 16)

set(gcf,'pos',[30 70 830 430])
% print('fig3','-depsc')

disp('done')

%=========================================================================
% figures 4, 11
%-------------------------------------------------------------------------
function figure4(ell,theta,sigma,pi,rho,eta,i,M,A,B)

figure()

subaxis(1,2,1, 'Spacing', 0.06, 'Padding', 0, 'Margin', 0.08);
eta = 0;
[PI,EC,LA,LB,Q0A,Q1A,Q0B,Q1Bn,ENn]=PI_mat(ell,theta,sigma,pi,rho,eta,i,M,A,B);
plot(PI,LA,'k-','Linewidth',3), hold on, plot(PI,LB,'k--','Linewidth',3)
axis square
ylim([0.85 1.55])
set(xlabel('$$\pi$$','interpreter','latex'), 'FontSize', 18)
set(gca,'FontSize',10.5,'TickLabelInterpreter','latex')
set(title('CRS','interpreter','latex'), 'FontSize', 15)

subaxis(1,2,2, 'Spacing', 0.06, 'Padding', 0, 'Margin', 0.08);
eta = 0.5;
[PI,EC,LA,LB,Q0A,Q1A,Q0B,Q1Bn,ENn]=PI_mat(ell,theta,sigma,pi,rho,eta,i,M,A,B);
plot(PI,LA,'k-','Linewidth',3), hold on, plot(PI,LB,'k--','Linewidth',3)
axis square
ylim([0.85 1.55])
set(xlabel('$$\pi$$','interpreter','latex'), 'FontSize', 18)
set(legend({'\,$$L_A$$','\,$$L_B$$'},'Interpreter','latex'), 'FontSize', 13)
set(gca,'FontSize',10.5,'TickLabelInterpreter','latex')
set(title('IRS ($$\eta=0.5$$)','interpreter','latex'), 'FontSize', 15)

set(gcf,'pos',[30 70 830 430])
% print('fig4','-depsc')

disp('done')

function figure11(ell,theta,sigma,pi,rho,eta,i,M,A,B)

figure()

subaxis(1,2,1, 'Spacing', 0.06, 'Padding', 0, 'Margin', 0.08);
eta = 0;
[PI,DELTA_A,DELTA_B]=sub_figure11(ell,theta,sigma,pi,rho,eta,i,M,A,B);
plot(PI,DELTA_A,'k-','Linewidth',3), hold on, plot(PI,DELTA_B,'k--','Linewidth',3)
axis square
xlabel('\fontsize{18} \pi')
ylim([0.0075 0.0130])
ax = gca; ax.YAxis.Exponent = -3;
set(gca,'FontSize',10.5,'TickLabelInterpreter','latex')
hT = title('CRS','interpreter','latex'); set(hT, 'FontSize', 15)

subaxis(1,2,2, 'Spacing', 0.06, 'Padding', 0, 'Margin', 0.08);
eta = 0.5;
[PI,DELTA_A,DELTA_B]=sub_figure11(ell,theta,sigma,pi,rho,eta,i,M,A,B);
plot(PI,DELTA_A,'k-','Linewidth',3), hold on, plot(PI,DELTA_B,'k--','Linewidth',3)
axis square
xlabel('\fontsize{18} \pi')
lT=legend({'\,$$\Delta_A$$','\,$$\Delta_B$$'},'Interpreter','latex'); set(lT, 'FontSize', 13)
ylim([0.0075 0.0130])
ax = gca; ax.YAxis.Exponent = -3;
set(gca,'FontSize',10.5,'TickLabelInterpreter','latex')
hT = title('IRS ($$\eta=0.5$$)','interpreter','latex'); set(hT, 'FontSize', 15)

set(gcf,'pos',[30 70 830 430])
% print('fig11','-depsc')

disp('done')

%=========================================================================
% figures 5, 6, 7, 8
%-------------------------------------------------------------------------
function figure5(ell,theta,sigma,pi,rho,eta,i,M,A,B)

figure()

global Bmin Bgap Bmax 
Bmin = 0.01; Bgap = 0.005; Bmax = 0.3; 

subaxis(2,3,1, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0;
sub_figure5(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0$$'   ,'interpreter','latex'), 'FontSize', 15)

subaxis(2,3,2, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.07;
sub_figure5(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.07$$','interpreter','latex'), 'FontSize', 15)

subaxis(2,3,3, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.14;
sub_figure5(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.14$$','interpreter','latex'), 'FontSize', 15)
set(legend({'\,$$L_A$$','\,$$L_B$$'},'Interpreter','latex'), 'FontSize', 15)

subaxis(2,3,4, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.21;
sub_figure5(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.21$$','interpreter','latex'), 'FontSize', 15)
set(xlabel('$$S_B$$','interpreter','latex'), 'FontSize', 15)

subaxis(2,3,5, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.28;
sub_figure5(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.28$$','interpreter','latex'), 'FontSize', 15)
set(xlabel('$$S_B$$','interpreter','latex'), 'FontSize', 15)

subaxis(2,3,6, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.35;
sub_figure5(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.35$$','interpreter','latex'), 'FontSize', 15)
set(xlabel('$$S_B$$','interpreter','latex'), 'FontSize', 15)

set(gcf,'pos',[30 70 1000 650])
% print('fig5','-depsc')

disp('done')

function figure6(ell,theta,sigma,pi,rho,eta,i,M,A,B)

figure()

global Bmin Bgap Bmax 
Bmin = 0.01; Bgap = 0.005; Bmax = 0.3; 

subaxis(2,3,1, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0;
sub_figure6(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0$$'   ,'interpreter','latex'), 'FontSize', 15)

subaxis(2,3,2, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.07;
sub_figure6(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.07$$','interpreter','latex'), 'FontSize', 15)

subaxis(2,3,3, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.14;
sub_figure6(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.14$$','interpreter','latex'), 'FontSize', 15)
set(legend({'\,Sell-probability','\,Buy-probability'},'Interpreter','latex','Location','southeast'), 'FontSize', 14)

subaxis(2,3,4, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.21;
sub_figure6(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.21$$','interpreter','latex'), 'FontSize', 15)
set(xlabel('$$S_B$$','interpreter','latex'), 'FontSize', 15)

subaxis(2,3,5, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.28;
sub_figure6(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.28$$','interpreter','latex'), 'FontSize', 15)
set(xlabel('$$S_B$$','interpreter','latex'), 'FontSize', 15)

subaxis(2,3,6, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.35;
sub_figure6(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.35$$','interpreter','latex'), 'FontSize', 15)
set(xlabel('$$S_B$$','interpreter','latex'), 'FontSize', 15)

set(gcf,'pos',[30 70 1000 650])
% print('fig6','-depsc')

disp('done')

function figure7(ell,theta,sigma,pi,rho,eta,i,M,A,B)

figure()

global Bmin Bgap Bmax 
Bmin = 0.04; Bgap = 0.001; Bmax = 0.08; 

eta = 0.2;
sub_figure5(ell,theta,sigma,pi,rho,eta,i,M,A,B)
axis square
ylim([0.8 1.4])
set(legend({'\,$$L_A$$','\,$$L_B$$'},'Interpreter','latex'), 'FontSize', 13)
set(gca,'FontSize',10.5,'TickLabelInterpreter','latex')
set(xlabel('$$S_B$$','interpreter','latex'), 'FontSize', 15)
xticks([0.04 0.05 0.06 0.07 0.08])
xticklabels({'0.04','0.05 ($$=S_A$$)','0.06','0.07','0.08'})

set(gcf,'pos',[30 70 415 430])
% print('fig7','-depsc')

disp('done')

function figure8(ell,theta,sigma,pi,rho,eta,i,M,A,B)

figure()

global Amin Agap Amax 
Amin = 0.01; Amax = 0.3; 

subaxis(2,3,1, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0;    Agap = 0.005; 
sub_figure8(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0$$'   ,'interpreter','latex'), 'FontSize', 15)

subaxis(2,3,2, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.07; Agap = 0.004;
sub_figure8(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.07$$','interpreter','latex'), 'FontSize', 15)

subaxis(2,3,3, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.14; Agap = 0.007;
sub_figure8(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.14$$','interpreter','latex'), 'FontSize', 15)

subaxis(2,3,4, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.21; Agap = 0.005; 
sub_figure8(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.21$$','interpreter','latex'), 'FontSize', 15)
set(xlabel('$$S_A$$','interpreter','latex'), 'FontSize', 15)

subaxis(2,3,5, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.28; Agap = 0.004;
sub_figure8(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.28$$','interpreter','latex'), 'FontSize', 15)
set(xlabel('$$S_A$$','interpreter','latex'), 'FontSize', 15)

subaxis(2,3,6, 'Spacing', 0.07, 'Padding', 0, 'Margin', 0.08);
eta = 0.35; Agap = 0.005; 
sub_figure8(ell,theta,sigma,pi,rho,eta,i,M,A,B)
set(title('$$\eta=0.35$$','interpreter','latex'), 'FontSize', 15)
set(xlabel('$$S_A$$','interpreter','latex'), 'FontSize', 15)

set(gcf,'pos',[30 70 1000 650])
% print('fig8','-depsc')

disp('done')

%=========================================================================
% figure 10
%-------------------------------------------------------------------------
function sub_figure10(ell,theta,sigma,pi,rho,eta,i,M,A,B,grid1,grid2)

row = grid1;
col = grid2;
Q00 = zeros(col,row);
Q11 = zeros(col,row);
[ELL,THETA] = meshgrid(linspace(0.01,0.99,row),linspace(0.01,0.99,col));
for j=1:numel(ELL)
    ell   = ELL(j);
    theta = THETA(j);
    [q0A,q1A,~,~,~]=symmetricEquilibrium(0.5,ell,theta,sigma,pi,rho,eta,i,M,A,B);
    Q00(j) = q0A;
    Q11(j) = q1A;
    fprintf('computing %.0f out of %.0f\n\n',j,numel(ELL))
end

AA = 2*(-1+ELL).*ELL.*THETA.*(1./Q00-1./Q11)./(-THETA+(-1+THETA)./Q11);
BB = -(-1+ELL).*THETA.*(-THETA+(-1+THETA)./Q00).*(-1./Q11.^2)./(THETA-(-1+THETA)./Q11).^2;
CC = -(1+(THETA-ELL.*THETA)./(-THETA+(-1+THETA)./Q11)).*(-1./Q00.^2);
JJ = (-1+1./Q11)./(Q00-Q11-log(Q00)+log(Q11));
KK = (-1+1./Q00)./(Q00-Q11-log(Q00)+log(Q11));
MM = 4*(-1+ELL).*ELL.^2.*THETA.*(-log(Q00)+log(Q11)+(Q00-Q11)./Q11)./(-THETA+(-1+THETA)./Q11);
NN = -(-1+ELL).*ELL.*THETA.*(Q00.*THETA-Q11.*THETA+log(Q00)-THETA.*log(Q00)+(-1+THETA).*log(Q11)).*(-1./Q11.^2)./(THETA-(-1+THETA)./Q11).^2;

NUMER = CC.*JJ.*MM + BB.*KK.*MM - 2*AA.*KK.*NN;

figure(), contour(ELL,THETA,sign(NUMER),[0,0],'k')
axis square, xticks([0.01 0.99]), yticks([0.01 0.99]), xticklabels({'0','1'}), yticklabels({'0','1'})
set(gca,'FontSize',20,'TickLabelInterpreter','latex')
set(xlabel('$$\ell$$','interpreter','latex'), 'FontSize', 20, 'position', [0.48 -0.02]) 
set(ylabel('$$\theta$$','interpreter','latex'), 'FontSize', 20, 'Rotation', 0, 'position', [-0.04 0.45]) 

function figure10(ell,theta,sigma,pi,rho,eta,i,M,A,B)

A = 0.03; B = A;
sub_figure10(ell,theta,sigma,pi,rho,eta,i,M,A,B,200,200)
set(gcf,'pos',[30 70 415 430])
text(.2,.835,sprintf('A'), 'Interpreter', 'latex', 'fontsize', 20);
text(.47,.47,  sprintf('B'), 'Interpreter', 'latex', 'fontsize', 20);
% print('fig10a','-depsc')

A = 0.06; B = A;
sub_figure10(ell,theta,sigma,pi,rho,eta,i,M,A,B,201,200)
set(gcf,'pos',[480 70 415 430])
text(.2,.835,sprintf('A'), 'Interpreter', 'latex', 'fontsize', 20);
text(.39,.39,  sprintf('B'), 'Interpreter', 'latex', 'fontsize', 20);
% print('fig10b','-depsc')

A = 0.1; B = A;
sub_figure10(ell,theta,sigma,pi,rho,eta,i,M,A,B,200,200)
set(gcf,'pos',[930 70 415 430])
text(.67,.67,sprintf('A'), 'Interpreter', 'latex', 'fontsize', 20);
text(.24,.24,  sprintf('B'), 'Interpreter', 'latex', 'fontsize', 20);
% print('fig10c','-depsc')

disp('done')
